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Recent breakthroughs in the experimental manipulation of strongly interacting atomic Rydberg 
gases in lattice potentials have opened a new avenue for the study of many-body phenomena. 
Considerable efforts are currently being undertaken to achieve clean experimental settings that 
show a minimal amount of noise and disorder and are close to zero temperature. A complementary 
direction investigates the interplay between coherent and dissipative processes. Recent experiments 
have revealed a first glimpse into the emergence of a rich non-equilibrium behavior stemming from 
the competition of laser excitation, strong interactions and radiative decay of Rydberg atoms. The 
aim of the present theoretical work is to show that local incoherent loss and gain of atoms can in 
fact be the source of interesting out-of-equilibrium dynamics. This perspective opens new paths for 
the exploration of non-equilibrium critical phenomena and, more generally, phase transitions, some 
of which so far have been rather difficult to study. To demonstrate the richness of the encountered 
dynamical behavior we consider here three examples. The first two feature local atom loss and gain 
together with an incoherent excitation of Rydberg states. In this setting either a continuous or a 
discontinuous phase transition emerges with the former being reminiscent of genuine non-equilibrium 
transitions of stochastic processes with multiple absorbing states. The third example considers the 
regime of coherent laser excitation. Here the many-body dynamics is dominated by an equilibrium 
transition of the “model A” universality class. 

PACS numbers: 32.80.Ee, 64.70.qj, 42.50.-p 


I. INTRODUCTION 

Identifying, classifying, and understanding the emer¬ 
gence of collective phenomena and other many-body ef¬ 
fects is a central objective of physics. In the past decades 
the refinement of experimental techniques for preparing, 
addressing and measuring atomic ensembles [T] opened 
entirely new possibilities for investigating not only sta¬ 
tionary, but also dynamical properties of quantum many- 
body systems. Rich non-equilibrium physics often stems 
from the presence of competing dynamical processes and 
strong interactions. Amidst different platforms, gases of 
atoms excited to high-lying Rydberg states are currently 
receiving increasing attention ms] as they feature con¬ 
siderable interactions via dipole-dipole or van der Waals 
forces. Their interplay with the laser excitation pro¬ 
cess — inducing coherent Rabi oscillations on the atomic 
populations — has shown to be the source of intricate 
behaviours, such as the formation of crystalline ground 
state structures CM], collectively-enhanced Rabi oscilla¬ 
tions coHia or the emergence of correlated equilibrium 
states UMI]. 

Since very recently, there has been considerable inter¬ 
est in understanding the many-body physics of interact¬ 
ing ensembles of Rydberg atoms in the presence of noise. 
While previously perceived as a detrimental feature, dis¬ 
sipation — caused by fluctuating atomic level shifts or 
radiative decay — can be in fact a source of an intrigu- 
ingly rich dynamics. Examples include the occurence of 
slow or glassy dynamics [TSHU], the relaxation into sta¬ 
tionary states with spatial correlations [22H24) , the obser¬ 
vation of intermittency and bistabilities [2H427] and the 


emergence of equilibrium EZHSI] and out-of-equilibrium 
universal behavior [32] . 

In this work we introduce a new scenario for the study 
of out-of-equilibrium phases and phase transitions with 
Rydberg atoms. The setting we have in mind consists of 
a background gas — acting as a large reservoir — from 
which Rydberg states are only excited at given spatial 
positions which are arranged in a regular lattice. Atoms 
from the background dynamically enter and leave these 
excitation spots. In conjunction with the laser-excitation 
and the strong interatomic interactions this local loss 
and gain dynamics leads to the emergence of non-trivial 
many-body dynamics. 

Such a scenario could be implemented in two rather 
different settings: in the first, a lattice of optical traps is 
immersed in a cold cloud of atoms and the traps are con¬ 
tinuously filled and depleted [33]. Current experiments 
aim at progressively slowing down this local dynamics 
by, for example, reducing the pressure of the background 
gas or increasing the strength of the optical confinement 
[341436] . These attempts could be reversed and, in prin¬ 
ciple, the setup could be “worsened” to the point that 
the timescale of the loss/gain dynamics becomes compa¬ 
rable with the other relevant dynamical processes. The 
second experimental setting consists of hot atomic gases 
confined in thermal vapour cells. Recently it has been 
shown that they indeed allow the observation of corre¬ 
lated many-body dynamics [20l [27| when Rydberg states 
are excited. Our envisioned setup is then realized by 
restricting the laser excitation to a regular array of ad¬ 
dressed spots. Thermal motion would move atoms in 
and out of these laser-illuminated regions, yielding the 
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desired loss/gain dynamics. 

Beyond introducing an additional dynamical process 
the consideration of local atom loss and gain might actu¬ 
ally relax a number of challenges that are currently faced 
by experimentalists when studying collective many-body 
behavior in dissipative Rydberg lattices. It might also 
simplify the modelling of Rydberg gases in which typi¬ 
cally radiative decay is accounted for as a dominant de¬ 
coherence mechanism: 


(i) It is not necessary to have (uniformly) determinis¬ 
tically loaded lattices, equal lattice confinement of 
ground state and Rydberg atoms or very low tem¬ 
perature states. 

(ii) It is not necessary to construct lattice potentials 
that equally trap ground state and Rydberg atoms. 
In fact it is not necessary to keep atoms trapped 
over an entire experimental run. 

(hi) One can employ very strongly interacting and high- 
lying Rydberg states that are typically long-lived. 
For such states the corresponding decay rate might 
simply be too small. In other words, it might be 
difficult to reach a regime in which the decay dy¬ 
namics is able to properly compete with the laser 
excitation and interatomic interactions, which thus 
almost entirely characterize the evolution. 

(iv) Even when acting on timescales that set it in com¬ 
petition with the driving, radiative decay is in¬ 
evitably accompanied by momentum kicks from 
photon recoil. Even when a Rydberg atom even¬ 
tually decays to the desired electronic ground the 
resultant heating might lead to loss of the atoms 
which can be accounted for in our description. 


Eor the sake of simplicity and in order to focus on the 
new aspects introduced by the loss/gain dynamics we will 
not consider radiative decay processes in this work. The 
underlying assumption is that the loss/gain dynamics is 
faster than that of the decay and/or that decay effectively 
induces a loss process via the mechanism described in 
point (v) above. 

The paper is organized as follows: after defining the 
model in Sec. [H] we introduce tw o sce narios that respec¬ 
tively feature a continuous [Sec. Ill and discontinuous 
[Sec. IV dynamical phase transition to an effective ab¬ 
sorbing subspace of states. Here we focus on a setting 
where the laser excitation of Rydberg states is described 
by a classical rate equation due to the presence of 
strong dephasing noise. In Sec. |V| we discuss the case 
of a coherent laser in the framework of a mean-field ap¬ 
proach. Here we show that the dynamics is described by 
the so-called ’Model A’ universality class, similar to what 
has recently been found in dissipative Rydberg gases with 
radiative decay [23 [29l [31] . Concluding remarks are pro¬ 
vided in Sec. IVIl 


II. THE MODEL 

We employ the standard description of a Rydberg lat¬ 
tice gas where each atom is modeled in terms of an effec¬ 
tive two-level system. The ground state |^) is coupled to a 
Rydberg nS'-state |t) through a laser with Rabi frequency 
ft and a detuning A with respect to the atomic transition 
(see Eig.j^a) for a visual representation). Within the ro¬ 
tating wave approximation the many-body Hamiltonian 
is then given by 

H = fl'^al +A'^nk + ^'^VkpUkUp, (1) 

k k k^p 

where Vkp = Cq/ \ Yk — represents the van der Waals 
(vdW) potential between pairs of excited atoms at posi¬ 
tions Yk and Fp, and the sum runs over all lattice sites 
k = 1... L. Interactions among ground-state atoms or 
between ground-state and Rydberg atoms are signifi¬ 
cantly weaker and will therefore be neglected. The op¬ 
erators refer to the Pauli matrices in the |t), |i) 

subspace, i.e. 

= Itfc) {ik\ + life) (tfel I 
<^k — {ik\ +i\ik) itkl , (2) 

^fe = life) (tfel - life) (tfel • 

Eurthermore, the local density of excitations is defined 
as n/c = |t/c) (t/c| and the density of ground state atoms 
as pk = \ik) 

In order to account for atom gain/loss in the lattice 
sites we add an effective third state |0), denoting an 
empty (or “inactive”) site. We also introduce the cor¬ 
responding local densities of active sites Ck = rik + Pk = 
|t/c) (t/c| + lik) This local loss and gain takes place 
with atoms from a background gas which is assumed to 
act as a bath. In other words, the surrounding cloud 
contains a much higher number of atoms than can be 
accommodated in the lattice and the recapture of a lost 
one is an unlikely event. Eirst of all, this suppresses 
correlations between loss and gain processes and allows 
us to treat them as being independent. Secondly, since 
atoms are constantly exchanged with new ones no corre¬ 
lations can be produced in the system via these processes. 
Thirdly, their occurrence probabilities are not apprecia¬ 
bly affected by the history of occupation of a given site, 
and can thus be considered Markovian. 

The relevant processes are schematically displayed in 
Fig.E] and summarized below: 

It) ^ |0) , It) ^ |0) , |0) ^ It). (3) 

with 7Dt, 7 d^ and 7 b^ being the corresponding rates. 
The first two processes describe the loss of a Rydberg 
and ground state atom, respectively. The third process 
corresponds to the capture of a ground state atom from 
the background gas. Note, that we do not consider the 
eventuality of Rydberg atoms being captured. The rea¬ 
son is that laser excitation to Rydberg states is restricted 
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(a) Internal Structure 



FIG. 1. Schematic representation of the system: (a,b) an optical lattice is realized within a background cloud of atoms. 
Atoms within the sites undergo laser-induced coherent transitions between their ground state W) and a high-lying (Rydberg) 
state It). The corresponding Rabi frequency and laser detuning are Q and A, respectively. A third state, |0), describes an 
empty or “inactive” site. Atoms are captured in and released from the sites with rates 7 b^ (capturing a ground-state atom), 
7Di (losing a ground-state atom) and 7Dt (losing an excited atom). The atomic states are furthermore subject to dephasing at 
a rate T and radiative decay from the Rydberg state into their ground state at rate k. (c) Rydberg atoms interact with a van 
der Waals potential Vkq, whose value for nearest neighbors is denoted by Vnn. The corresponding energy shift of the Rydberg 
state in the vicinity of an excited atom is sketched, (d) Emergent branching processes in the limit of strong dephasing: in the 
case analysed in Sec. |III| (square lattice) a single excitati on h as the potential for branching (with rate A) which enables the 
production of larger clusters. In the case discussed in Sec. IV (triangular lattice) the system parameters are chosen such that 
two nearby excitations are required for a cluster to grow. 


to local sites and consequently Rydberg atoms are not 
produced in the background gas. Hence, the transition 
|0) ^ It) could only occur if a Rydberg atom is cap¬ 
tured which had been previously expelled from another 
site, which is unlikely. Note, that we are also neglect¬ 
ing processes which lead to the occupation of a given 
site with multiple atoms. In fact, in experiments with 
microtraps such multi-occupancies are suppressed due to 
the collisional blockade [38] |39] . In circumstances where 
such suppression is not taking place the so-called dipole 
or Rydberg blockade [Mllinilll] is ensuring that each site 
can only feature a single Rydberg excitation. This also 
limits the local site dynamics to a restricted space at the 
expense of a possibly varying local (density-dependent) 
Rabi frequency [42]. For the sake of simplicity we will 
not consider such situations in this work. 

In addition to the loss/gain dynamics we consider the 
presence of noise, which dephases local superpositions 
between the states |t) and ||) at a rate F. The origin of 
this noise can be fluctuating background fields that result 
in random atomic level shifts, the broadening of atomic 
lines due to Doppler broadening m or interaction effects 
[43] . or a spectrally broad excitation laser [2]. 

In the presence of the described coherent and dissipa¬ 
tive processes the evolution of the density matrix p of 
the system is governed by a quantum master equation in 
Lindblad form [44] [45] 


dtp = -i [H, p] + '^ji 

i,k 


Li,kpL\f. 


1 

2 






^/eP^/c 


1 

2 


{nk,p} 


(4) 


Here the index i runs over the three different sources of 


noise introduced in Eq. Q, while k over all lattice sites; 
{A, B} = AB + BA is shorthand for an anticommutator. 

= |0/e) U/c|, = \^k) (t/e| and = \ik) (0/c | 

are the corresponding local jump operators; the last term 
accounts for dephasing, implemented via the jump oper¬ 
ators n/c- 


III. CONTINUOUS TRANSITION FROM/TO 
AN ABSORBING STATE 


We start by considering a situation in which the de¬ 
phasing amplitude F is much larger than the Rabi fre¬ 
quency and the other dissipative rates (F ^ 

In this regime, the dynamics is effectively described by 
means of a classical stochastic equation [I8]|37l|46]|47|: 
the underlying separation of timescales permits the adi¬ 
abatic elimination of the portion of the phase space sub¬ 
ject to dephasing. Correspondingly, the evolution of the 
density matrix p of the system is projected onto the 
dissipation-free subspace [484150] . which in this case cor¬ 
responds to the sole diagonal components in the basis 
(i.e. a basis of classical spin configurations) [37]. At the 
leading order in a perturbative expansion in powers of 
7 i/F and U/F the truncated density matrix p evolves 
according to 


dtp = 'B,^k {(^kP^k - (^kp) + 

k 

4 " ^i,kpk.,k ~ 2 


( 5 ) 
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where 


A. Mean-field approach 


Afe = 



q^k 


n„ 


-1 


(6) 


is a configuration-dependent rate. Defining the diago¬ 
nal part of the Hamiltonian {H Iq^q) “classical” 

component, the second term in the brackets of Eq. (§ 
corresponds to the square of the classical energy change 
accompanying a spin-flip at site k. Spin-flips that result 
in a significant increase or decrease in energy are there¬ 
fore strongly suppressed. In the scenario investigated in 
this section we choose the detuning A such that it is 
opposite to the interaction energy Enn between neigh¬ 
boring excited atoms (A = — Enn) [see Fig. [^c)]. Hence, 
exciting an atom right next to an isolated already ex¬ 
cited one incurs no energy difference and therefore oc¬ 
curs at the maximal rate = X = 4D^/r (see 

Fig. 0 (d)). We further assume to be working in a regime 
where —A = ^ F, i.e. the interaction surpasses the 

dephasing strength. In this regime any atom that has 
more than 1 excitation in its neighborhood (or none at 
all), can only change its internal state at a rate of order 
A/c oc D^r/A^. This rate is significantly smaller than 
A^max) such proccsscs are strongly suppressed. 

For brevity, we shall refer to all of them as “off-resonant” 
processes. Among them, we emphasize that the creation 
of excitations in an empty neighborhood, occurring at a 
rate Ak ~ D^F/A^, is also strongly suppressed. 

In order to gain some first insight into the expected 
many-body dynamics we assume for a moment that all 
off-resonant processes can be neglected. In this regime 
Eqs. and §) display features that characterize a class 
of stochastic processes m, such as the contact pro¬ 
cess and other branching-annihilating ones, which are 
known to undergo genuine non-equilibrium phase transi¬ 
tions [52l |53] . These systems feature an absorbing phase 
with strictly zero density n = {uk) /T, i.e. in the ab¬ 
sence of excitations no new ones can be created. Further¬ 
more, the production of excitations can only proceed via 
nucleation of clusters and is, in this sense, local (it cannot 
occur arbitrarily far from pre-existing excitations). It is 
important to note that in the present case the absorb¬ 
ing phase does not consist of a unique state, but rather 
an absorbing space that is spanned by the entire set of 
configurations of sites which are either in state 0 or state 
In general there is dynamics taking place within the 
absorbing manifold as all the absorbing configurations 
can be visited via the interplay of the local loss and gain 
processes with rates 7 d^ and 7 b^. This requires us to 
consider also the density of active sites r] = (e/^) jL 

when analyzing the dynamics of the system. 


The mean-field approximation discards correlations 
between different sites — i.e., for every local observable 
O/c we substitute (OkOp) (Ok) {Op) if k ^ p — and 
permits the formulation of closed equations of motion 
for the expectation values of the densities of excitations 
n and of active sites 7 : 

dtn = + E ( i 

i=o ^ ^ 


dtV = 7Bi — (7 b^ + 7d^) 7 + (Td^ — 7Dt) 


(7a) 

(7b) 


Here, 2 ; is the lattice coordination number (number of 
nearest neighbours per site) and Xj is shorthand for the 
rate of a flipping process occurring in the presence of 
j excited neighbors. This means that Ai characterizes 
the resonant processes introduced above (and Ai = A = 
4D^/F), whereas the remaining values refer to the off- 
resonant processes and read 


- 








A. ( 8 ) 


In a first approximation we neglect them altogether and 
Eq. (7a) becomes 


dtfi = —7Dt^ + (7 “ 2n) (1 — n) 


z-l 


(9) 


This equation together with Eq. ( |7b| ) predicts a transi- 
tion from the region A < Ac = 7Dt(7B4. + 
which admits only the absorbing solution n = 0 to the 
region A > Ac in which the system displays a finite den¬ 
sity n > 0 in the long-time limit. In Fig. we report the 
corresponding phase diagram for the choice A = — 64F, 
7Dt = 7Di = 7 d in the X-f] plane, where the symbol 
7 denotes the stationary density of active sites. The 
threshold value Ac identifies the critical point of a con¬ 
tinuous transition between the two phases. To its right, 
the density scales linearly (n ^ A — Ac), while at the 
critical point its value decays to 0 in time according to 
the power-law n{t) ^ Ijt. The density of active sites re¬ 
laxes to the finite value 7 b^/( 7 b^ + 7 d^)- Consequently, 
at the mean-field level, this system undergoes a transi¬ 
tion which shows some of the characteristic features of 
directed-percolation (DP) universality [5T] . 

Let us now discuss the role of the off-resonant terms. 


Those with j > 1 in Eq. (7a) do not affect the fundamen¬ 
tal properties of the transition, as they vanish for n ^ 0. 
Therefore, they can only shift the position of the critical 
point according to the relative statistical weights Xj . The 
j = 0 term, on the other hand, constitutes a relevant — 
albeit small — perturbation that brings the system away 
from the critical point. The reason is that it accounts 
for production of excitations in an empty neighborhood 
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FIG. 2 . Stationary density of ex cita tions n extracted from 
the mean-field equations and ( [7b| for —A = Vnn = 64r, 
7 b^ = 0.0IF and yot = 7 d^- The data is shown as a function 
of the branching rate A for resonant processes and the den¬ 
sity of active sites f] = 7Bi/(7B^ + 7Di)- The color scale is 
bounded by Umax = 0.5. The red dashed line corresponds to 
the values taken by the critical rate Ac for different values of 
ff. A cross section is displayed in the inset for 77 = 0.8 (along 
the cyan horizontal line in the main figure), which highlights 
the mean-field scaling behavior n ~ A —Ac. The green, dashed 
line corresponds to the same curve calculated including the 
leading off-resonant processes relevant in a Rydberg gas . As 
expected, the introduction of the latter makes the transition 
smoother, but deviations are only visible in close vicinity to 
the critical point Ac. 


and thus prevents the aforementioned subspace of config¬ 
urations devoid of excitations from being strictly absorb¬ 
ing. This term smooths the transition into a crossover, 
as highlighted in the inset of Fig. The magnitude of 
this effect can be suppressed by increasing the detuning 
A. When sufficiently small it allows the observation of 
the mean-field scaling behavior for values of A > Ac- 


B. Numerical analysis 


In order to investigate the effect of fluctuations which 
are not captured by the mean-field treatment we per¬ 
form numerical Monte Carlo simulations, using a state 
in which all sites are occupied with a Rydberg atom as 
the initial condition. We set the rates 70 ^ = 7 Dt = 7d, 
iBi = O.OIF, Vnn = 64F = —A and collect data for 
several values of the parameters 7 d and Q. For this par¬ 
ticular choice of parameters the loss and gain dynam¬ 
ics decouples from the excitation dynamics. This can be 
seen directly in Eq. (7b) which is valid beyond mean-field. 


Consequently, the density of active sites 77 reaches expo¬ 
nentially fast (on a timescale ( 73 ^ + 7 d)~^) the steady- 
state value f) = 7 Bt/ ( 7 Bt + 7 d^ 

For Rydberg gases one needs to account for the fact 
that the off-resonant production of excitations and the 
long-range tails of the vdW potential affect the emer¬ 
gence of the phase transition. As we discuss further be¬ 


low, these features actually constitute a source of addi¬ 
tional noise which to some extent may obscure the an¬ 
ticipated scaling behaviors. In order to shed light on the 
fundamental critical properties of the transition we have 
therefore also simulated a dynamical process in which we 
replace the first term of the r.h.s. of Eq. ©by 


E 

k,ie{k} 


^ni{a+^ia^ - PkfJ-) , 


( 10 ) 


with {k} denoting the set of nearest-neighboring sites of 
site k and d being the dimension. After this replace¬ 
ment we have a pure branching process (as found, e.g., 
in the contact process mentioned above m) producing 
excitations from nearby ones at a rate A/(2d). The nor¬ 
malization by the quantity 2d — which corresponds to 
^ for square lattices — is meant to compensate for the 
fact that in this case multiple excitations enhance the 
rate. We emphasize that, although different, the two 
processes we consider share the same fundamental prop¬ 
erties: the absorbing subspace is the same and, apart 
from off-resonant events, branching is the only way to 
increase the number of excitations. Eurthermore, in the 
presence of low-densities — as happens in the proximity 
of the critical point — the action of the branching terms 
in Eqs. ^ and (10) is analogous up to multiplicative fac¬ 
tors. Eor brevity, in the following we shall refer to the 
new stochastic process as the “pure” instance and to the 
original process [Eq. <§] as the “Rydberg” one. 

In Eig. I^we show the stationary value (obtained from 
Monte Carlo simulations) of the density of excitations n 
in the fj-X plane for ID and 2 D square lattices in both 
cases. The pure and the Rydberg processes display quali¬ 
tatively the same behavior. As expected, in the pure case 
the transition from the absorbing to the active phase is 
significantly sharper. Beyond that, two interesting fea¬ 
tures can be observed. Eirst, the simulations seem to 
suggest that the critical point Ac diverges as the station¬ 
ary density of active sites 77 decreases and that below 
a certain threshold fjc the transition disappears entirely. 
Second, there is a qualitative difference in the static and 
dynamic scaling behavior when varying 77 . The stationary 
properties remain unaffected and always display, within 
numerical accuracy, a scaling behavior n r\j (A —Ac)^ with 
a critical exponent p compatible with the DP one for 
both one- and two-dimensional processes (see Eig. 1^. In 
contrast, the dynamical approach to stationarity changes 
continuously. This means that when approaching the 
threshold value 77 c, the critical exponent of the algebraic 
decay n(t) ^ smoothly decreases from a value which, 
in ID, is comparable with the one of pure DP to 0. 

The latter feature is strongly reminiscent of the behav¬ 
ior of stochastic processes with multiple absorbing states 
as reported in Refs. imii], which provide a qualitative 
explanation of our observations. Even though in our case 
the absorbing space is not made of individually-absorbing 
states, the excitation dynamics effectively perceives them 
as such, since it stops completely as soon as the first ab¬ 
sorbing configuration is reached. Moreover, in the cases 
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FIG. 3. Phase diagrams of the pure and the Rydberg process (see text) in the ry-A plane for a ID chain of 100 sites and a 2D 
square lattice of 20 x 20 sites. The parameters are chosen as —A = Rnn = 64r, 7b^ = O.OIF and 7d^ = 7Dt = 7 d. The color 
scale is set with respect to the maximal value the density can take, i.e., nmax = 1 for the pure process and 1/2 for the Rydberg 
one. Numerically computed exponents /3 (static) and S (dynamic) are displayed in the panels. The selected parameter ranges 
are shown as a cyan and a green line on the main plot. For the ID pure process we also show (in log-log scale) the critical 
profiles of the stationary density n as a function of Aa = A — Ac (lower-left inset) and of its evolution in time (upper-right 
inset) to highlight the scaling behavior. For comparison we provide the known DP exponents [51]: /3 id = 0.276, /32d = 0.584, 
Sin = 0.159, and fco = 0.451. 


discussed in Refs. iniEi the dynamic exponent is also 
not constant but instead varies continuously as a func¬ 
tion of the initial conditions, e.g. the initial density. In 
our simulations we start from a fixed initial condition (all 
atoms present and excited). However, the fast loss/gain 
dynamics rapidly constructs an “effective initial condi¬ 
tion” with an active site density ry determined by the 
rates 7 b^ and 70 - Since this initial condition is varied 
under a change of 7 b^ and 7 d this might be a possible 
explanation for the observed variation of the dynamic 
exponent. 

The Rydberg case features off-resonant processes and 
thus displays a smoothed transition. Moreover, it re¬ 
quires stronger driving for the active phase to appear. 
This can be understood by noting how clusters of excita¬ 
tions actually hinder their own growth. For instance, if 
we consider a pair of nearby excitations, elongating it to 
a three-excitation segment (tti^ttt) faces the presence 
of next-nearest-neighbor interactions. Because of them, 
the rate at which this process occurs is no longer A but 
instead given by 


a(NNN) 






. 1/2 
'^NNN 


( 11 ) 


Our choice of parameters, Rnnn = f^NN/2^ = F, implies 
= A/5 and hence the branching rate is effectively 
reduced. Further growth along the same direction ex¬ 
periences much smaller corrections and thus continues 
approximately at a rate A/5. The situation worsens if 
we consider branching orthogonally with respect to the 
two original excitations, since in this case the distance 
between next-nearest neighbors is reduced to \/2 times 
the lattice spacing, implying Rnnn = f^NN/(A/2)^ = 8F 
and yielding an effective rate = A/257. This ex¬ 

plains the suppression of the stationary density in the 
2D case with respect to the ID one. The relevance of 
this effect can be drastically reduced by partially remov¬ 
ing the tails of the vdW interactions using a microwave 
dressing scheme (see Refs. [32l [55H57] ). In other words, 
by coupling two Rydberg levels with a strong microwave 
field one can obtain a hybridization of the relative in¬ 
teractions. For appropriate choices, the latter features 
a crossover threshold separating a short-distance regime 
displaying the usual vdW decay from a long-distance one 
which is instead suppressed with respect to the previ¬ 
ous one and can be considered approximately flat. This 
length-scale thus acts as a cutoff for the potential tails. 
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FIG. 4. Stationary density of excitations n (we select the 
stable solution with highest value) extracted from the mean- 
field equations ( |13[ ) and ( |7b[ ) . The data is shown as a function 
of A (branching rate) and f] (stationary density of active sites) 
for -A = 2yNN = 64r, 7b^ = O.OIF and yot = Here 
^max = 0.5. The dashed line corresponds to the critical rate 
Ac for different values of r/. In the inset we plot the profiles 
of the three stationary solutions of Eq. (13) for fj ~ 0.65 in a 
parameter interval indicated by the cyan horizontal line in the 
main figure. The absorbing solution is displayed as a solid, 
cyan line which never leaves 0. The upper solid orange line 
represents the stable solution with finite density which only 
occurs for A > Ac. The dashed green line corresponds to the 
third, unstable solution. 


IV. DISCONTINUOUS ABSORBING-STATE 
PHASE TRANSITION 

The continuous phase transition discussed above is not 
merely due to the competition between system-filling and 
system-emptying processes. It also strongly relies on the 
resonance condition requiring the presence of a single ex¬ 
citation nearby, which constitutes a fundamental connec¬ 
tion with other classical branching processes m- 

To demonstrate this we consider a 2D triangular lat¬ 
tice {z = 6) of Rydberg atoms which is irradiated by 
a laser with detuning — A = 2 Unn- This implies that 
atoms now require exactly two excited neighbours to be 
effectively brought in resonance with the laser [see Fig.[^ 
(d)]. With this constraint, the geometrical structure of 
the lattice becomes relevant: for instance, clusters can¬ 
not grow on a square lattice, since neighboring sites do 
not share common neighbors. 


A. Mean-field approach 


The density of active sites is evolving according to 
Eq. (7b) in this case as well. The evolution equation 


for the excitation density n can be cast in the form (7a), 
but with Xj redefined as A 2 = A and 


- 






{^y + A^j-2f A^j-2y 


< A. (12) 


Hence, the dominant contribution is now given by the 
j = 2 term, i.e.. 


dtn « -jD^n + 2 ^ 

(13) 

Introducing the parameter x = 2z{z — 1)Q ‘^= 
Xz{z — l)/27Dt, these mean-field equations predict a dis¬ 
continuous transition from a phase x < Xc in which the 
only acceptable stationary solution is n = 0 to an ac¬ 
tive phase (x > Xc) in which two further solutions ap¬ 
pear (constituting a saddle-node bifureation [58], see inset 
of Fig. 1^. These correspond to additional real roots of 
the r.h.s. of Eq. ( [T^ and are identified by the equality 
~ 2n-) = (1 — . Thus, the threshold value for 

X corresponds to 

Xc = niin — 73 — ^ ' -r—3. (14) 

ne[o,i/2] n{r] — 2n) (1 — nY~^ 

Branching off from a common point (the value n yield¬ 
ing the minimum in Eq. ([l4|)), one of these solutions is 
unstable (dashed green line in the inset of Fig. un¬ 
der small perturbations and decreases asymptotically to 
0, whereas the other one (solid blue line) is stable and 
increases up to 7^/2. The absorbing solution remains al¬ 
ways stable. A phase diagram for this case is shown in the 
main panel of Fig. where we display, for various choices 
of A and r/, the largest stable mean-field solution, in or¬ 
der to highlight the finite jump experienced by this value 
when the boundary Ac = 27 DtXc/(^(^ — 1)) is crossed. 
The quantitative effect of the off-resonant terms is again 
barely noticeable at the mean-field level and we thus do 
not display it. We recall, however, that the subspace of 
configurations with n = 0 is not perfectly absorbing and 
thus the lower branch is slightly lifted from 0. 


B. Numerical analysis 

To study the system beyond mean-field we employ 
Monte Carlo simulations. The parameters are identi¬ 
cal to the previous simulations with the exception that 
f^NN = — A/2 = 32F. Again, we set 7d^ = Tot = 7d 
and use a completely filled initial state from which the 
fast site dynamics will generate an effectively random 
initial configuration. In Fig. we display the station¬ 
ary density n in the (2-7 d plane. Note, that for the 
sake of clarity we do not include interactions beyond 
nearest neighbours, i.e. we assume that the tails of the 
vdW potential are modified for example by the previ¬ 
ously mentioned microwave dressing jSailSHST]. Their 
inclusion does not qualitatively change the phase dia¬ 
gram, although it makes it more difficult to highlight the 
discontinuous nature of the transition. As predicted by 
the mean-field study we observe the presence of both an 
absorbing phase and an active one. At a first glance the 
data suggest the presence of a continuous crossover in¬ 
stead of a discontinuous jump. This can be reconciled 





















n1.0 


- 0.5 


0.0 


FIG. 5. (a) Stationary density profile in the X-ij plane for the branching process defined on a 25 x 25 triangular lattice. The 
color scale extends up to Umax = 0.5. As in the mean-field treatment in Fig. we fix —A = 2yNN = 64F, = O.OIF and 

7Dt = 7 d^ = 7d. For a sufficiently large population of trapped atoms (sufficiently small loss rate 7 d) and large enough A 
(oc the process becomes clearly capable of maintaining a finite density of excitations for long times, (b) Counting statistics 
p{n) of the final density n(tmax) at fixed Q = 0.125F and Ftmax = 25000 as a function of fj. Bimodality is a signature of the 
discontinuity of the phase transition and becomes clearly visible in the inset section taken along a section at ry ~ 0.7194. 


with the previous Subsection’s conclusions by consider¬ 
ing that in the presence of two stable stationary solutions, 
different realisations of the stochastic process end up in 
either one or the other. However, this information is lost 
once the average is taken. 


We have therefore computed the counting statistics — 
or probability distribution — of the excitation density n 
at the maximal time of our simulations, Ft^ax = 15000. 
It is displayed in Fig.[^b) and clearly highlights the dis¬ 
continuous first-order nature of the transition showing 
that indeed, for a certain range of parameters, the dis¬ 
tribution becomes bimodal. This confirms that some ini¬ 
tial conditions decay into empty configurations, whereas 
other ones maintain a finite density n > 0 for long times. 


V. “MODEL A” UNIVERSALITY 


For completeness we now briefly discuss the situation 
in which dephasing is not strong enough in order to war¬ 
rant a description of the dynamics of the system in terms 
of effectively classical evolution equations. In this situ¬ 
ation the full quantum master equation 0 needs to be 
solved which can be done numerically only for systems 
consisting of a few sites. To nevertheless gain a quali¬ 
tative understanding of the emerging phase structure we 
employ a mean-field treatment. The single-site expecta¬ 
tion values (n/c) = n, {pk) = p, {cf%) = and {cf\) = 


evolve according to equations 

fi = — 7Dt^ (15a) 

P =- 7DiP +7B4.(1 --P) (15b) 

Sx = _(A + Vn)Sy _ ^ + (15c) 


Sy = 2n{p - n) + (A + Vn)S^ - ^ + ^Dt + ^Dl gy^ 

(15d) 

where V = 14,(1 - 5kq). 

First, we focus on investigating the stationary proper¬ 
ties. Setting the time derivatives to zero one can relate all 
steady-state expectation values to the excitation density 
n, i.e. 

sy = (16a) 


P = 


-7Dt«’+ 7 b 4,(1 - n) 


= -2 


7B4. + 7D4. 
A + Vn 


7Dt 


r + 7Dt + 7 d4, 


n. 


(16b) 


(16c) 


This leads a closed polynomial equation for the excitation 
density in the stationary state which reads 

n [a + 6(c + n)^] = 1, (17) 

where we have introduced the parameters 

_ 27B4. + 7Dt + 7D4, (r + 7Dt + 7D4.)7Dt(7B4, + 7D4.) 


7B4. 


4H^7B4. 


fvV 7Dt(7Bi + 7D4.) 

V H / 7b4. (r + 7Dt + 7D4.) 


and c 


A 

V' 


( 18 ) 
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FIG. 6. Sta tionary p hase diagram of the mean-field equations 
of motion ( 15a|15d] ) for three different detunings A in the 
plane spanned hy Q./V and = 7Dt/^ = 7 d/V. The 

shaded areas correspond to the bistable regions for F = 0.0ly 
and 7Bi = 7 d- Points lying outside represent choices for 
which there is only one physical solution. From top to bot¬ 
tom, the detuning is hxed at —0.2551/, —0.22V and —0.15V. 
The boundaries of each shaded region meet with vanishing 
net angle at the corresponding critical point (red disc). The 
solid black line displays the path taken by the critical point 
as A is varied. The dashed lines show how the prohle of 
this path shrinks when the dephasing rate F is increased. 
From the outermost to the innermost, the values correspond 
to F = O.OOiy, O.Oiy (solid line), 0.04y and O.OSV. The 
transition disappears when F/V > 7Bi/4(7B^ + 7d) = 0.125. 
The bistable region is also absent when A/V lies outside the 
interval [-97Bi/[16(7Bi + 7 d)], 0] = [-9/32,0]. 


Similar mean-field equations have been obtained in a 
number of recent works m EH EH im EU] which inves¬ 
tigate the dynamics of Rydberg gases in the presence of 
radiative decay. A known feature of Eq. is that it de¬ 
scribes a bistable behavior, not fundamentally different 
from the one encountered in the previous Section, i.e. for 
certain values of the physical parameters the mean-field 
equations of motion admit two, instead of one, stable 
stationary solutions. 

Figure [^provides a representative example of the phase 
diagram’s structure in the 7 d“^ plane, obtained in the 
plot for the specific choice 70 ^ = 7Dt = 7Bt = 7 d- 
The bistable regions are shown as shaded areas. In gen¬ 
eral, the lines delimiting these domains meet with van¬ 
ishing angle in a critical point (see also the discussion 
in Ref. [28]) which for fixed parameter c is located at 
Uc = — 9/8c, be = — 27/8c^. This identifies a unique di¬ 
rection, which in the parameter space spanned by a, 5 and 
c reads Sa = -Sbe^jO where Sa = a — ac and Sb = b — be- 
Varying the parameters in such a way to follow this path, 
the profile of the stationary density displays a branching 
at the critical point which can be related to the sponta¬ 
neous breaking of a Z 2 symmetry [28j. The solutions for 


the stationary density of excitations take here the form 
n = He ^ Sn^ with the critical density ric = — 2c/3 and 


6n = 


0 _ {6b < 0) 


(19) 


Crossing the critical point along any other direction 
yields instead a behavior 6n ^ 6b^^^, which relates to 
an explicit breaking of the symmetry. This is identical to 
the situation that is encountered in the thermodynamics 
of a ferromagnet in the presence of a finite magnetic field. 

To understand the dynamical behavior near the criti¬ 
cal point we first note that for our choice of parameters 
7Dt = 7Dt = 7 d the site dynamics decouples. The corre¬ 
sponding equation is 


dt{n+p)=r] = ^Bi-{lBi+7D)v- ( 20 ) 


Linearising the remaining equations in the proximity of 
the critical point, we find two attractive modes and a van¬ 
ishing one, signalling the presence of a critical slowing- 
down of the dynamics, i.e., an algebraic relaxation to¬ 
wards the stationary value n{t) — ric ^ t~^. 

Recent experiments in Rydberg gases suggest that a 
mean-field treatment such as the one presented here can 
correctly capture qualitative dynamical and static fea¬ 
tures [saisi]. A recently-developed variational method 
EiEn], which improves on the mean-field, showed in this 
context that the bistable region is in fact ’replaced’ by a 
first-order-transition line terminating in a critical point. 
The latter is related to an emerging universality belong¬ 
ing to the so-called “model A” (or Ising-Glauber) class 
j6T] . The lines separating the bistable and stable regions 
can consequently be interpreted as spinodal lines, re¬ 
sulting in the appearance of long-lived metastable states 
which have been observed experimentally m and whose 
properties relate to the mean-field predictions [28] . 

Hence, the results in this section indicate that the in¬ 
terplay between the coherent laser excitation, the inter¬ 
action and the loss/gain dynamics can drive the system 
towards a classical equilibrium critical point. Finally, we 
remark that here — contrary to the cases presented in 
Secs. |III| and[TV]— the presence of dephasing is inessential 
for the emergence of criticality. Moreover, the transition 
disappears when the dephasing rate exceeds 


Fc 


1 _TBt_ 

2 27b^ + 7Dt + 7Dt 


V 


- 7Dt - 7Dt 


as indicated in Fig. From a physical point of view, 
this can be understood by the fact that the various cases 
reported here rely on facilitation effects induced by the 
blue-detuning of the laser with respect to the atomic 
transition. As discussed above, this means that the pres¬ 
ence of excited atoms enhances the production of new 
ones (at a certain distance). This happens due to the 
fact that the interactions effectively shift the atomic lev¬ 
els. Therefore, if the broadening of the latter induced by 
dissipation becomes comparable with the shift itself (e.g., 
F V), the facilitation effect disappears and so does the 
critical behavior. 
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VI. CONCLUSIONS 

In this work we have illustrated a number of non¬ 
equilibrium phenomena that can be explored within lat¬ 
tice gases of Rydberg atoms in the presence of local loss 
and gain processes. As a matter of fact, the latter, which 
are currently regarded as unwanted sources of noise to 
overcome, turned out to be key for the emergence of the 
collective dynamics reported above. 

In the limit of strong dephasing the Rydberg gas real¬ 
izes several instances of transitions from and to absorbing 
states. Facilitating the flipping of atoms which possess 
a single excited neighbor yields a continuous transition 
whose dynamic properties vary smoothly depending on 
the average density of active sites, while the static ones 
remain unchanged. This is strongly reminiscent of the 
behavior of stochastic processes with multiple absorbing 
states. On the other hand, on a triangular lattice and 
with a resonance condition asking for two neighbors to 
be excited the system undergoes a discontinuous phase 
transition, highlighted by the bimodal structure of the 
stationary density of excitations. 

In the limit where the coherent laser excitation dom¬ 
inates the dephasing processes the static and dynamical 
properties of the system are determined by an equilib¬ 
rium critical point belonging to the “model A” univer¬ 
sality class. This hints at the possibility of investigating 
these collective behaviors in a setting complementary to 
the ones in which these phenomena were previously stud¬ 


ied, all of which relied on the presence of radiative decay. 

A particularly appealing aspect of the present study is 
the prospect to explore non-equilibrium phase transitions 
with absorbing states. Although the underlying univer¬ 
sality classes - all of which relate, to a certain degree, 
to directed percolation — have been thoroughly investi¬ 
gated in the past, there are not many condensed matter 
systems which are known to display the corresponding 
critical behaviors [62f[64] . Rydberg gases in the presence 
of local loss and gain permit the detailed exploration of 
directed percolation universality in all three dimensions 
j32] and also in the presence of a manifold of absorbing 
states. Peculiar features such as the dependence of the 
critical dynamics on the initial state, which have been 
studied in the framework of idealized model systems, can 
now be in principle explored experimentally. 
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